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EHD ponderomotive forces and aerodynamic flow control using plasma actuators 
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We present a self-consistent two-dimensional fluid model of the temporal and spatial development 
of the One Atmosphere Uniform Glow Discharge Plasma (OAUGDP®). Continuity equations for 
electrically charged species Ng , Nj, O^, O^" and electrons are solved coupled to the Poisson equa- 
tion, subject to appropriate boundary conditions. It was used an algorithm proposed by Patankar. 
The transport parameters and rate coefficients for electrons at atmospheric pressure are obtained 
by solving the homogeneous Boltzmann equation for electrons under the hydrodynamic assumption. 
Operational variables are obtained as a function of time: electric current; surface charge accumu- 
lated on the dielectric surface; the memory voltage and the gas voltage controlling the discharge. It 
is also obtained the spatial distribution of the electric field, the populations of charges species, the 
resulting ponderomotive forces, and the gas speed. 

PACS numbers: 47.65.+a,52.80.Pi,47.70.Fw,47.70.Nd,51.50.+v,47.62.+q 

Keywords: Electrohydrodynamics, Dielectric barrier discharges (DBD), atmospheric pressure glow discharge 
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I. INTRODUCTION 

The development of the One Atmosphere Uniform 
Glow Discharge Plasma (OAUGDP®) has made it possi- 
ble to generate purely electrohydrodynamic (EHD) pon- 
deromotive (body) forces 0, 0, 0, 3 . Such forces are 
generated without a magnetic field and with small in- 
tensity currents crossing the plasma. In fact, only RF 
displacement currents produce the body forces that ac- 
celerate the plasma. Two methods were devised for flow 
acceleration 0- E| : 1) Peristaltic flow acceleration and 
2) Paraelectric flow acceleration. Only the last method 
is analyzed in this work. Paraelectric flow acceleration 
is the electrostatic analog of paramagnetism: a plasma 
is accelerated toward increasing electric field gradients, 
while dragging the neutral gas with it. Applications span 
from propulsion and control systems in aeronautics, to 
killing any kind of bacterium and virus (see Ref. pj). 

The role of plasma in aerodynamic research has been 
increasing, since it constitutes a significant energy mul- 
tiplier modifying the local sound speed and thus lead- 
ing to modifications of the flow and pressure distribution 
around the vehicle 0, H, • Plasma actuators have been 
shown to act on the airflow properties at velocities below 
50 m/s [13 ■ 

In default of a complete model of a OAUGDP® reac- 
tor, Chen ^3 built a specific electrical circuit model 
for a parallel-plate and coplanar reactor, modeling it 
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as a voltage-controlled current source that is switched 
on when the applied voltage across the gap exceeds the 
breakdown voltage. 

Although there is still lacking a detailed characteriza- 
tion of such plasma actuators, with only boundary layer 
velocity profiles measured using a Pitot tube located 1-2 
mm above the flat panel |12| being available, we present 
in this paper a self-consistent two-dimensional modeling 
of temporal and spatial development of the OAUGDP® 
in an "airlike" gas. 



II. NUMERICAL MODEL 
A. Assumptions of the model 

We intend to describe here the glow discharge regime, 
with emphasis on flow control applications of the plasma. 
Gadri [l3| has shown that an atmospheric glow dis- 
charge is characterized by the same phenomenology as 
low-pressure dc glow discharges. 

No detailed plasma chemistry with neutral heavy 
species is presently available; only the kinetics involving 
electrically charged species supposedly playing a deter- 
minant role at atmospheric pressure: Nj, , Oj, O^, 
and electrons, is addressed. The electronic rate coeffi- 
cients and transport parameters are obtained by solving 
the homogeneous electron Boltzmann equation under the 
hydrodynamic regime assumption When obtaining 
the charged species populations as well the electric field 
controlling their dynamics, the following electrohydrody- 
namics (EHD) effects are studied in this work: 

• ponderomotive forces acting on the plasma hori- 
zontally and perpendicularly to the energized elec- 
trode; 

• gas velocity. 
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FIG. 1: Schematic representation of electrode geometry of an 
energized OAUGDP plasma panel. 

The simulation domain is a 2-DIM Cartesian geometry 
(see Fig. 2J with total length along the X-axis Lx = 0.4 
cm and height Ly = 0.4 cm; the width of the dielectric 
surface along the X-axis is 0.3 cm in Case Study I and 
0.1 cm in Case Study II. The dielectric relative permittiv- 
ity, supposed to be a ceramic material, is assumed to be 
e r = 100; the dielectric thickness is in all cases set 0.065 
cm. The capacity of the reactor is determined through 
the conventional formula Cd s — £o£ r S/d. The electrodes 
thickness is supposed to be negligible. 

B. Transport parameters and rate coefficients 

The working gas is a " air like" mixture of a fixed frac- 
tion of nitrogen (Sn 2 — [N 2 ]/N = .78) and oxygen 



{$o 2 — [0 2 ]/N — 0.22), as is normally present at sea 
level at p = 1 atm. 

The electron homogeneous Boltzmann equation |14j is 
solved with the 2-term expansion in spherical harmonics 
for a mixture of -/V2 — 22% O2. The gas temperature is 
assumed constant both spatially and in time, T g = 300 
K, and as well the vibrational temperature of nitrogen 
T V (N 2 ) = 2000 K and oxygen T v (0 2 ) = 2000 K. The 
set of cross sections of excitation by electron impact was 
taken from [l5j . 

At atmospheric pressure the local equilibrium assump- 
tion holds: Transport coefficients (v^, ^fon^ Mei Mpj D e , 
D p ) depend on space and time (r, t) only through the lo- 
cal value of the electric field E(r, t). This is the so called 
hydrodynamic regime. 

Ion diffusion and mobility coefficients were taken from 
[HI, [i -.N = 6.85 x 10 21 V" 1 m" 1 s" 1 (on the range of 

E/N with interest here), fi +N = 6.91 x 10 21 V" 1 m" 1 

s-\ and fi N +N = 5.37 x 10 21 V" 1 m" 1 s" 1 . 

The reactions included in the present kinetic model are 
listed in Table HJ It is assumed that all volume ioniza- 
tion is due to electron-impact ionization from the ground 
state and the kinetic set consists basically in ionization, 
attachment and recombination processes. The kinetics 
of excited states and heavy neutral species is not consid- 
ered. 

To obtain a faster numerical solution of the present 
hydrodynamic problem it is assumed that the gas flow 
does not alter the plasma characteristics and is much 
smaller than the charged particle drift velocity. This as- 
sumption allows a simplified description of the flow. For 
more concise notation, we put n p2 = [N^ ]; n p \ = [N^]] 
n p = [0£]; n„ = [O2 ], and n e = [e]. 

The balance equations for N^" (at atmospheric pressure 
the nitrogen ion predominant is N4") is: 

—^-+V-{n p2 v p2 ) = S N2 N 2 n pl K lcl -(3Nn p2 -K r2 n p2 n e . 

(1) 

The balance equation for is: 



—jjj- + V • (n p iv p2 ) = n e vf ^ + K ic2 [N 2 ]n p2 - (3un n n p i - f3n e n p i - K icl [N 2 ] 2 n pl . (2) 

I 

The oxygen ion considered is O 2 and its resultant balance introduced and its balance equation was written as: 
equation is given by 

Q n -Kf + v ' K v «) = v att n e - Pun pl n n - K d n p n n . (4) 

+ V • (n p v p ) = n e v io 2 n - j3 %i n n n p - (3n e n p . (3) 

Finally, the balance equation for electrons can be written 
As oxygen is an attachment gas, the negative ion was in the form: 
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dn e , . 

-Qj- + V • (n e v e ) = n e [y. 



N 2 



+ v ?on ~ v att) ~ 0n e (n p + rtpi) + K d [0 2 ]n n - K r2 n pl n e 



(5) 



r 



To close the above system of equations we use the drift- 
diffusion approximation for the charged particle mean 
velocities appearing in the continuity equations: 



(6) 



where fii and Di represent the charged particle mobil- 
ity and the respective diffusion coefficient. The applied 
voltage has a sinusoidal wave form 



V(t) = V dc + V ahx(wt), 



(7) 



where Vdc is the dc bias voltage (although here we fixed to 
ground, Vdc = 0) and u is the applied angular frequency. 
Vq is the maximum amplitude with the root mean square 
voltage in this case of study V rms = 5 kV and the applied 
frequency / = 5 kHz. 



The total current (convective plus displacement cur- 
rent) was determined using the following equation given 
by Sato and Murray 



h{t) = 



V 



On? 

p dz 



D, 



dn e 



D, 



dn n 



F, L dv 



V 



f dE L 
v \ dt 



E L dv, 



(8) 



where J v dv is the volume occupied by the discharge, 
is the space-charge free component of the electric field. 
The last integral when applied to our geometry gives the 
displacement current component 



dv. 



(9) 



Auger electrons are assumed to be produced by impact 
of positive ions on the cathode with an efficiency 7 = 
5 x f0~ 2 , so that the flux density of secondary electrons 
out of the cathode is given by 



Jse{t) = 7Jp(*)i 



(10) 



with j p denoting the flux density of positive ions. In 
fact, this mechanism is of fundamental importance on 
the working of the OAUGDP®. 

Due to the accumulation of electric charges over the di- 
electric surface, a kind of "memory voltage" is developed, 
whose expression is given by: 



Vrn(t) = ^- / I d (t')dt' 
<^ds Jtn 



+ V m (t ) 



(11) 



Here, Cd s is the equivalent capacitance of the discharge. 

The space-charge electric field was obtained by solving 
the Poisson equation 



AV = (n p 

£0 



n„). 



(12) 



The boundary conditions are the following: 



electrode (Dirichlet boundary condition): V(x,y 
0,t) = V-V m ; 



• dielectric (Neumann boundary condition) 



The flux of electric charges impinging on the dielectric 
surface builds up a surface charge density a which was 
calculated by balancing the flux to the dielectric 



da 
Iff 



— = e(ir„ 



|r e ,„|). 



(13) 



Here, r p „ and r e)ra represent the normal component of 
the flux of positive and negative ions and electrons to the 
dielectric surface. Furthermore, it is assumed that ions 
and electrons recombine instantaneously on the perfectly 
absorbing surface. 

The entire set of equations are solved together, at each 
time step, self-consistently. 



III. METHOD OF RESOLUTION OF FLUID 
EQUATIONS 

The particle's governing equations are of convection- 
diffusion type. They are solved usin g a method pro- 
posed by Patankar [ljj (see also Ref. [2(|). According 
to this method, let L(<p, d(f), d 2 cj), ...) = S be a homo- 
geneous differential equation in 0, with a source term 
S. Then the procedure of Patankar consists in replacing 
L{(j),d^),d 2 (j), ...) = S by apipp = ^2 k a k (j) k + b, where P 
is the central point of the mesh. 

The chosen time step is limited by the value of the di- 
electric relaxation time. For the present calculations the 
total number of computational meshes used is (100x100). 
This fair condition allows calculating an entire cycle with 
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TABLE I: List of reactions taken into account in our model. Rate coefficients were taken from Ref. [r?|. 



kind of reaction 


Process 


Ionization 




Op _L /V"+ 


Ionization 


e + 2 -> 


2e + 0+ 


3-body electron attachment 


e + 2 +0 2 


-> + o 2 


3-body electron attachment 


e + 2 + N 2 


-> o 2 ~ + iV 2 


Collisional detachment 


Oz +o 2 ~ 


-> e + 20 2 


e-ion dissociative recombination 


N+ + e 


-» 2iV 


e-ion dissociative recombination 


0+ + e 


-» 20 


2-body ion-ion recombination 


2 + N+ - 


O2 + iV 2 


Ion-conversion 


N+ + N 2 + N 2 -> N+ + e 


Recombination 


N+ + e 


-> 2AT 2 


Ion-conversion 


N+ + N 2 ^ 


N+ + 2N 2 



Rate coefficient 



ion 
,Pl a 



Kai = 1.4 x 10- 29 (f2)exp(-600/T 9 )A'i(r 9 ,re) (cm 6 /s) 
K a2 = 1.07 x lb- 31 (^l) 2 ^2(T 9 ,r e ) (cm 6 /s) c 

K d = 2.7 x 10- 10 Y^exp(-5590/r 9 ) (cm 3 /s) 

= 2.8 x 10~ 7 y|^ (cm 3 /s) 

/? = 2.8 x l(r 7 ./f& (cm 3 /s) 



2 x 10-\/f^|l + 10 



-«iV(fO) 2 ] (cm 3 /s) 



Kiel = 5 x 10" 29 (cm 6 /s) 
if r2 = 2.3 x 10~7(Te/300) ' 56 
K ic2 = 2.1 x 10~ 16 exp(T 9 /121) (cm 3 /s) 



"Data obtained by solving the quasi-stationary, homogeneous elec- 
tron Boltzmann equation. See Ref. 1 1 -J for details. 
'With K x = exp(700(T e - T s )/(T e T s )) 
c With K 2 = exp(-70/T 9 ) exp(1500(T c - T g )/{T e T g )) 



an Intel Pentium 4 (2.66 GHz) in a reasonable CPU time 
of about 30 hours per cycle, limiting to a reasonable value 
the relative error \&w e /w e \, with w e designating the elec- 
tron drift velocity. Stationarity was attained typically 
after 4-5 cycles. Equations ^ EH are integrated succes- 
sively in time supposing the electric field is constant dur- 
ing each time step, obtaining a new value of the electric 
field after the end of each time step. The method used 
to integrate the continuity equations and Poisson equa- 
tion was assured to be numerically stable, constraining 
the time step width to the well known Courant-Levy- 
Friedrich stability criterion. 



IV. RESULTS 

The simulations were done for a two-dimensional flat 
staggered geometry, as sketched in Fig. ^ This is essen- 
tially a 'surface discharge' arrangement with asymmetric 
electrodes. It is assumed that the plasma is homogeneous 
along the OZ axis. 



Electrical characteristics 




FIG. 2: Electric current, applied voltage, gas voltage and 
memory voltage as a function of time. Conditions: Case I. 



Solid curve: current; dot curve: 
dashed curve: V. 



V m \ dashed-dot curve: Vo 



In Fig. [21 it is given the evolution along a period of 
time of the calculated electric current, applied voltage, 
gas voltage and memory voltage. The OAUGDP®, and 
as well generically a DBD, occurs in configurations char- 
acterized by a dielectric layer between conducting elec- 
trodes. At about 740 Volts, electron avalanches develop, 
replenishing the volume above the surface with charged 
particles. Hence, the charged particles are flowing to the 
dielectric (see Ea. ll3f) start accumulating on the surface, 
and build-up an electric field that prevents the occur- 
rence of a high current, and quenches the discharge de- 



velopment at an early stage. 



B. Electrical field and potential 

Fig-Elshows the electric field during the first half-cycle 
at the instant of time t — 1.9 x 10 -5 s. The energized elec- 
trode is the anode and the electric field follows Aston's 
law, its magnitude remaining on the order of 10 5 V/cm at 
a height of 8 x 10~ 5 m above the electrode and attaining 
lower magnitude above the dielectric surface, typically 
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FIG. 4: Calculated space averaged ponderomotive forces per 
unit volume as a function of time. Solid curve: F x ; dot curve: 
F y . Case study I 



FIG. 3: Electric field along OX and OY near the energized 
electrode at time t = 1.9 x 10 s at first half-cycle. Condi- 
tions: Case study I, with V = 5 kV, / = 5 kHz. 

on the order of 10 3 V/cm. The electric field magnitude 
is strongest in region around the inner edges of the en- 
ergized electrode and dielectric surface (which is playing 
during this half-cycle the role of a pseudo-cathode). 

During the avalanche development a strong ion sheath 
appears. In fact, as the avalanche develops an ion sheath 
expands along the dielectric surface until it reaches the 
boundary. With the ion sheath travels an electric field 
wave, with some similarities with a solitary wave. The 
speed of its propagation in the conditions of Fig. is 
about 150 m/s. See Refs. 0, 21] for very elucidating 
explanation of this phenomena. 



C. Paraelectric gas flow control 

The theory of paraelectric gas flow control was devel- 
oped by Roth . The electrostatic ponderomotive force 
Fe (units N/m 3 ) acting on a plasma with a net charge 
density p (units C/m 3 ) is given by the electrostatic pon- 
deromotive force and can be expressed under the form 



F E = -e V£ 2 . (14) 

In order to verify whether electrostriction effects could 
be playing any significant role, it was also calculated the 



electrostriction ponderomotive force 

Fes = --e E 2 Ve r . (15) 

Here, e r — 1 — vi is the relative permittivity of the 
plasma, v en is the electron-neutral momentum transfer 
frequency and u p is the plasma frequency. We found 
that this force term is negligible, contributing at maxi- 
mum with 1 % to the total ponderomotive force. Sub- 
sequently, the ponderomotive forces were averaged over 
the area of calculation. Comparing the calculated space 
averaged ponderomotive forces per unit volume shown in 
Figs. BJ El it is seen that when the electrode width in- 
creases they become one order of magnitude higher. On 
average, during the second half-cycle the ponderomotive 
force magnitude decreases. This happens when the volt- 
age polarity is reversed and the energized electrode play 
the role of cathode. This is due to a reduction of the 
potential gradient on the edge of the expanding plasma 
(see also Ref. HJ). Calculations of EHD ponderomotive 
force have shown that its maximum intensity is attained 
during electron avalanches, with typical values on the or- 
der of 5 x 10 9 N/m 3 . F x points along OX (propelling 
direction), while F y points downwards (boundary layer 
control). 

D. Gas speed 

Using Bernoulli law (see Ref. |]|) it can be obtained 
the induced neutral gas speed 

v = E,f^= J-F X L X . (16) 

V p Vp 
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clear the slight decrease of the gas speed during the cath- 
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FIG. 5: Calculated space averaged ponderomotive forces per 

unit volume as a function of time. Solid curve: F x ; dot curve: FIG. 7: Space averaged gas speed as a function of time. Case 
F y . Case study II II 
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FIG. 6: Space averaged gas speed as a function of time. The 
straight solid line is a linear fit showing the increase of gas 
speed with time. Case I 

Here, F x is the calculated space average ponderomotive 
forces per unit volume, and p = 1.293 Kg/m 3 . Fig. 
shows the gas speed along the entire cycle in Case I. The 
average value of the gas speed is around 15 m/s while the 
experimental value, measured with a Pitot tube 1-2 mm 
above the surface, is 5 m/s as high as is for nearly the 
same operational conditions |l2^ . 

As can be seen in Fig. the gas speed increases to 
about 20 m/s when the dielectric surface decreases. It is 



ode cycle. This is related to the decrease of ponderomo- 
tive forces, as discussed above. As we assumed charged 
particles are totally absorbed on the dielectric surface, 
the swarm of ions propagating along the dielectric sur- 
face are progressively depleted, dwindle with time. How- 
ever, it is worth to mentioning (see also Ref. [2^) that 
in certain conditions the inverse phenomena can happen, 
a bigger dielectric width feeding up the ion swarm with 
newborn ions and thus inducing an increase of the gas 
speed. How long its width can be increased is a matter 
of further study. 



V. CONCLUSION 

A 2-DIM self-consistent kinetic model has been imple- 
mented to describe the electrical and kinetic properties of 
the OAUGDP®. It was confirmed that the electric field 
follows the Aston's law above the energized electrode. 
EHD ponderomotive forces on the order of 5 x 10 9 N/m 3 
can be generated locally during the electron avalanches, 
their intensity decreasing afterwards to values well below 
on the order of 10 4 -j- 10 5 N/m 3 . On the cathode side the 
EHD ponderomotive forces can decrease 1.5-5-2 orders of 
magnitude, due probably to a smaller important poten- 
tial gradient. The ponderomotive forces (and as well the 
gas speed) tend to increase whenever the energized elec- 
trode width augments relatively to the dielectric width. 

This code will help to design an advanced propulsion 
system, achieving flow control in boundary layers and 
over airfoils by EHD means, with numerous advantages 
over conventional systems. 



7 



[1] J. R. Roth, Physics of Plasmas 2003 10 2117 
[2] C. Liu, J. R. Roth, Paper 1P-26, Proceedings of the 
21st IEEE International Conference on Plasma Science, 
Santa Fe, NM, June 6-8 1 994, ISBN 7803-2006-9, pp. 
97-98 

[3] J. R. Roth, D. M. Sherman, and S. P. Wilkinson, AIAA 
Journal 2000 38 1166 

[4] J. R. Roth, P. P.-Y. Tsai, C. Liu, M. Laroussi, and P. 
D. Spence, "One Atmosphere Uniform Glow Discharge 
Plasma", U. S. Patent # 5414324, Issue May 9 (1995) 

[5] Roth, J R, Industrial Plasma Engineering, Vol.1: Princi- 
ples (Bristol,IOP,1995) 

[6] Roth, J R, Industrial Plasma Engineering, Vol. 2: 
Application to Nonthermal Plasma Processing (Bris- 
tol,IOP,2001) 

[7] P. Bletzinger, B. N. Ganguly, D. Van Wie and A. 

Garscadden, J. Phys. D: Appl. Phys. 2005 38 R33-R57 
[8] Alfredo Soldati, Sanjoy Banerjee, Phys. Fluids 1998 10 

1742 

[9] W. Shyy, B. Jayaraman, and A. Andersson, J. Appl. 

Phys. 2002 92 6434 
[10] Jerome Pons, Eric Moreau and Gerard Touchard, J. 

Phys. D: Appl. Phys. 2005 38 3635 
[11] Zhiyu Chen, IEEE Trans. Plasma Sci. 2003 31 511 



[12] J. Reece Roth, R. C. M. Mohan, Manish Yadav, Jozef 
Rahel and, Stephen P. Wilkinson, AIAA Paper 2004-0845 

[13] R. B. Gadri, IEEE Trans. Plasma Set. 1999 27 36 

[14] C. M. Ferreira, L. L. Alves, M. Pinheiro, and P. A. Sa, 
IEEE Trans. Plasma Sci. 1991 19 229 

[15] Siglo Data Base: http:cpat.ups-tlse.fr 

[16] R. S. Sigmond, Gas Discharge Data Sets for Dry Oxy- 
gen and Air, Electron and Ion Physics Research Group 
Report, the Norwegian Institute of Technology, The Uni- 
versity of Trondheim (1979) 

[17] Kossyi I A, Kostinsky A Yu, Matveyev A A, and Silakov 
V P, Plasma Sources Sci. Technol. 1992 1 207 

[18] R. Morrow and N. Sato, J. Phys. D: Appl. Phys. 1999 
32 L20-L22 

[19] S. V. Patankar, Numerical heat transfer and fluid flow, 

(New York: Taylor & Francis, 1980) 
[20] N. Pinhao, Modeling of the Discharge in Halogen 

Quenched Geiger-Muller Detectors (in Portuguese), Ph.D 

Thesis, Technical University of Lisbon, 1997 
[21] J. P. Boeuf and L. C. Pitchford, J. Appl. Phys. 2005 97 

103307 

[22] A. Shvydky, V. P. Nagorny and V. N. Khudik, J. Phys. 
D: Appl. Phys. 2004 37 2996 



